First principles calculation of vibrational Raman spectra in large systems: signature 

of small rings in crystalline Si02 
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We present an approach for the efficient calculation of vibrational Raman intensities in periodic 
systems within density functional theory. The Raman intensities are computed from the second 
order derivative of the electronic density matrix with respect to a uniform electric field. In 
contrast to previous approaches, the computational effort required by our method for the evalu- 
ation of the intensities is negligible compared to that required for the calculation of vibrational 
frequencies. As a first application, we study the signature of 3- and 4-membered rings in the 
the Raman spectra of several polymorphs of Si02 , including a zeolite having 102 atoms per unit cell. 



PACS numbers: 71.15.-m,78.30.-j,71.15.Mb 

Vibrational Raman spectroscopy [Q is one of the most 
widely used optical techniques in materials science. It 
is a standard method for quality control in production 
lines. R is very effective in determining the occurrence 
of new phases or structural changes at extreme condi- 
tions (high pressure and temperature), where it is often 
preferred to the more difficult and less readily available x- 
ray diffraction experiments based on synchrotron sources 
Moreover, it can be used in the absence of long- 
range structural order as for liquid or amorphous mate- 
rials |§, p. The theoretical determination of Raman 
spectra is highly desirable, since it can be used to asso- 
ciate Raman lines to specific microscopic structures. 

Density functional theory (DFT) |^ can be used to 
determine with high accuracy both frequencies and in- 
tensities of Raman spectra. Vibrational frequencies can 
be efficiently determined nsvag first order response ^. 
Within this approach Raman intensities (RI) calculation 
is also possible, but requires a computational time signifi- 
cantly larger and is not practical for large systems. Thus, 
while many examples of frequency calculations have been 
reported so far , RI were predicted from first-principles 
in a very limited number of cases involving systems with 
a small number of atoms j9[ |l^, |ll|. In this Letter we 
show that it is possible to obtain RI in extended solids 
with a computational cost negligible with respect to that 
required for the frequency determination. The efficiency 
of our approach will lead ab-initio calculations to become 
a routine instrument for the interpretation of experimen- 
tal Raman data. Our method is based on second order 
response to DFT. In particular, we compute the second 
order derivative of the electronic density matrix with re- 
spect to a uniform electric field, using pseudopotentials 
and periodic boundary conditions. As a first application 
we calculate Raman spectra of several Si02 polymorphs, 
including a zeolite having 102 atoms per unit cell jl^ . 

In a Raman spectrum the peak positions are fixed 
by the frequencies Uu of the optical phonons with null 
wavevector. In non-resonant Stokes Raman spectra of 



harmonic solids, the peak intensities I'^ can be computed 
within the Placzek approximation p| as: 
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r oc |e, ■A'^-e,!^ — (n, + l). 
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where (gs) is the polarization of the incident (scat- 
tered) radiation, Ui, = {exp{huJi,/kBT) — 1)^^, T is the 
temperature, and 



(2) 



Here f ' is the electronic energy of the system, Ei is the 
Z*'' Cartesian component of a uniform electric field, u/c-y 
is the displacement of the 7*'* atom in the k^^ direction, 
Af^ is the atomic mass, and vf^^ is the orthonormal vi- 
brational eigenmode v. 

Linear response [Q, ^ can be used to determine lOu, 
v^r", and also the dielectric tensor *e°° defined as = 
5irn — {'^T^ /^)d'^£^^ / {dEidEm)^ whcrc VI is the cell volume. 
RI have been computed |l^ through Eq. (|l|), obtain- 
ing A" by finite-differences derivation of *e °° with respect 
Ufc^. This approach requires 36 A'^^* linear response cal- 
culations, where A^^* is the number of atoms. Thus, the 
scaling of the RI calculation is the same as that of the 
frequency calculation with a much larger prefactor. This 
has limited the applications of this approach to small 
systems. RI have also been computed from the the dy- 
namical autocorrelation functions of 'e°° in a molecular 
dynamics (MD) run ||ll[. This approach also copes with 
liquids or anharmonic solids, but is very demanding, re- 
quiring the calculation of e*°° at each MD step. 

Alternatively, RI can be obtained knowing the sec- 
ond order derivative of the DFT density matrix p = 
X^i) li being jV't,) the normalized occupied Kohn- 

Sham (KS) eigenstates In fact, according to the well 
known Hellmann-Feynman theorem 

2 Tr <^ p 
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where Tr{0} is the trace of the operator O, and v*"^* is the 
external ionic potential (the KS self-consistent potential 
is V^^ = \/Hxc _|_ ycxt^ where jg ^j^g g^j-|^ -j-j^g 

Hartrec and the exchange-correlation potential). Thus 



= 2Tr 



dEidErnduk^ " y\dEidE„^J duk^ 

The p/{dEidEjn) calculation requires six second-order 
calculations, instead of the 36 A^^' first-order calculations 
needed for the finite-differentiation Because of this 
better size-scaling, the A"^ calculation through Eq. (|^) 
is much more efficient and the time for RI calculation is 
negligible compared to that for lUi, in large systems. 

The approach based on Eq. has already been used 
in isolated molecules but never in extended sys- 
tems. Indeed, in solids the calculation of d^p/ (OEidEm) 
is not trivial because the position operator, required by 
the electric field perturbation, is ill-defined in periodic 
boundary conditions. Because of this, although a for- 
malism to calculate derivatives of p at any order was pro- 
posed by Gonze already in 1995 only very recently 
Nunes and Gonze were able to include perturbations 
due to macroscopic electric fields. To do that, they use 
the polarization-Berry phase formalism , arguing that 
this concept remains valid in the presence of finite elec- 
tric fields. This approach has been applied so far to a one 
dimensional non-self-consistent model . In the follow- 
ing we give an expression for the second derivative of p, 
that does not require the Berry phase formalism to cope 
with uniform electric fields, and we use it to compute A'^ 
in real systems with a DFT self-consistent Hamiltonian. 

The derivative of p with respect to two generic pertur- 
bation parameters A and p is: 

d^p _ ^ 
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where P = (1 — p) is the projector on the empty state 
subspace, the sums over v and v' run over the occupied 
states, and Irii'^''^^) are the second derivatives of the occu- 
pied KS-orbitals in the parallel-transport gauge ||]. Ac- 
cording to our derivation: 



' dx ' 



Gy 



dV 



KS 



dX 

Q2yKS 

dXdp 



\i'v), 



dp ' dX 



dX ' dp 



(5) 



(6) 



Here, 
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TABLE I: Raman activity in Si computed with our approach 
(Tsor), and by finite differences (7fd). N is the number of 
inequivalent k-points. 
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is the Green function operator projected on the empty 
states \'ipc) [@, [A, B] = AB - BA, and the first deriva- 
tive of the density matrix is: 



dp 



(7) 



Since dV^^/dX and d^V^^/{dXdp) depend on dp/dX, 
dp/ dp, and d'^p/ {dXdp), Eqs. should be solved self- 

consistently. 

The advantage of the present formulation, compared 
to that of Ref. 1^, lies in the introduction of the com- 
mutators of Eqs. (HH). Thanks to the commutators, all 
the quantities needed with our formalism are well de- 
fined in an extended insulator, even if the perturbation 
p or A are the component i?; of a a uniform electric field, 
i.e. if dV^^/dX = -en + dV^'^'^/dEi (l|, being n 
the Z*'' Cartesian component of the position operator r, 
and e the electron charge. In particular, in an insulator, 
the commutators [r, p] and [r, dp/ dp] in Eqs. (||,||) are 
well defined, bounded operators, since the density ma- 
trix is localized ((r"|p|r') goes to zero exponentially for 
|r"-r'| ^oo). 

Finally, in a periodic system, the right-hand side of 
Eq. (||) can be easily computed in terms of the \v^), that 
are the periodic parts of the Bloch-wavefunctions 
with reciprocal-lattice vector k, using the substitutions: 



n, 



dp 



,k|5lO(w!J'i, ,k 



dki 
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dE, 



^' dh 
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where I and m are Cartesian indexes, c is an empty band 
index, v and v' are occupied band indexes, and Pk is the 
projector on the empty subspace of the point k. In our 
implementation, the derivative with respect to ki in the 
right-hand side of Eq. (||) is computed numerically by 
finite-differences, using an expression independent from 
the arbitrary wavefunction-phase, as in Refs. p^ . 

We test our approach on Si in the diamond phase, 
where the Raman activity is determined by 7 = adefl/du 
1^, where a = 10.20 a.u. is the lattice spacing and u the 
displacement of one atom along the (1,1,1) direction [ p^ . 
We compute 7 for various grids of k-points, using both 
our second order response method (7sor) ^md by finite 
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FIG. 1: Vibrational Raman spectra of various Si02 poly- 
morph powders. Measurements are from Refs. |^^. Theo- 
retical frequencies are rescaled by +5%, and the spectra are 
convoluted with a uniform Gaussian broadening having 4.0 
cm~^ width. 



differentiation with respect to the atomic displacement 
(Tfd), Tab. H At convergence the two approaches are 
completely equivalent. 

As a second application, we consider tetrahedral Si02- 
In this class of materials, that includes the all-silica zeo- 
lites, the quartz, cristobalite, tridymite and coesite poly- 
morphs of Si02, and vitreous sihca (v-Si02), each Si 
atom is tetrahedrally coordinated to four O atoms and 
each O atom is bonded to two Si atoms. The properties 
of these systems can be effectively described in terms of 



the n-membered rings (n-MRs) of tetrahedra contained 
in their structure ||- E.g., a clear correlation be- 

tween the presence of 3- and 4-MRs and the degradation 
of optical v-Si02 fibers under UV radiation has been ob- 
served Q. In the v-Si02 Raman spectra the two sharp 
peaks at 490 cm^^ {Di line) and 606 cm^^ {D2 line), 
have been attributed to the breathing mode (BM) of the 
O atoms towards the ring center of 4-MRs and 3-MRs, 
respectively This attribution has been confirmed by 
DFT vibrational frequency calculations . The attribu- 
tion would be further supported by experimental mea- 
surements on well characterized crystalline polymorphs 
containing 3- and 4-MRs. However, the strong Raman 
peak at 520 cm"-'^ in coesite, a phase that contains 4- 
MRs, is shifted by 30 cm~^ with respect to the Di line 
in v-Si02, and no Raman measurements has been pub- 
lished on the H-ZSM-IS zeolite, that is the only known 
Si02 crystalline polymorph with 3-MRs |l^. Interest- 
ingly this zeolite contains 4-MRs as well. 

To clarify this topic, we compute the Raman spectra of 
a-quartz, coesite, a-cristobalite, and H-ZSM-18 p0[ . 
In Fig. ^ we compare our results with the available exper- 
imental spectra . The vibrational frequencies are sys- 
tematically underestimated by 5% by our calculation. To 
simplify the comparison with the experiments, in Fig. |l| 
and|^, the theoretical frequencies are multiplied by a scal- 
ing factor of 1.05. The ability of the method in reproduc- 
ing quantitatively all the measured features is evident. 

In order to associate Raman peaks of Fig. |l| to the 
small-ring BMs, we project the vibrational eigenmode 
on the subspace generated by the BMs of a given kind 
of rings, TZ, and on the corresponding complementary 
subspace, 1Z. We use the two resulting projected vectors 
to decompose A'' so that A" — A^ + A^. Since /"^ is 

quadratic in A^ sec Eq. @), I" = H^J- + Ceriap, 

where /overlap term bihnear in A^ and A^. A 

Raman peak can be associated to a ring BM (i.e. the 
Raman activity is solely due to the BM ) if, and only if, 

^ I ^overlap I- 

The structure of H-ZSM-lS ||l^ contains two equiva- 
lent 3-MRs and two kinds of 4-MRs which we will call 
4-MRso, and 4-MRsi |||. In Fig. |, we show the pro- 
jected Raman spectra of the zeolite and the coesite. In 
the H-ZSM-18 spectrum, the peaks at 485, and 615 cm^^ 
are very well described by the BM of 4-MRso and 3-MRs, 
respectively. A direct analysis of the vibrational eigen- 
modes shows that both BMs are decoupled from other 
modes. The frequencies of the two peaks are very close 
to those of the measured Di and D2 lines in v-Si02 (490, 
and 606 cm~^), thus confirming that these lines are due 
to rings BMs ||^, |^. However, the presence of small-MRs 
in a structure does not guarantee, in general, the occur- 
rence of completely decoupled BMs. This is the case of 
the 4-MRs in coesite and the 4-MRsi in the zeolite, whose 
BMs exhibit a large |/ovcriapl' ^ig. |[ These over- 
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FIG. 2: Raman intensities projected on the breathing modes 
of various rings labeled 3-MRs, and 4-MRS3; (see the text). 
For clarity, the overlap intensity (/overlap in the text) is shifted 
vertically. 



laps imply the existence of a coupling with other modes, 
that, in turn, explains the 30 cm~^ difference between 
the 4-MRs frequency of coesite and that of the Di line of 
v-Si02- A comparable frequency shift from the Di line 
is observed, with opposite sign, for the 4-MRsi BMs in 
the zeolite. 

In conclusion, with the aim of building an instrument 
for the routine interpretation of Raman spectra, we de- 
veloped a method for the efficient calculation of Raman 
intensity. We computed the Raman spectra of Si02 poly- 
morphs containing up to 102 atoms. We found that: i) 
not all the small-membered-rings have decoupled breath- 
ing modes, ii) the H-ZSM-IS zeolite provides decoupled 
breathing mode of 4- and 3-membered rings, whose fre- 
quencies nicely coincide with the Di and D2 lines of vitre- 
ous silica. An experimental determination of the Raman 
spectra of this zeolite can thus provide an experimental 
calibration for the determination of the density of decou- 
pled small membered rings in vitreous silica. 
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